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ABSTRACT 

A high resolution finite element method for the solution of problems involving high 
speed compressible flows is presented. The method uses the concepts of flux-corrected 
transport and is presented in a form which is suitable for implementation on completely 
unstructured triangular or tetrahedral meshes. Transient and steady state examples 
are solved to illustrate the performance of the algorithm. 
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INTRODUCTION 


Over the past few years, there has been an ongoing interest in the application 
of unstructured grid finite element methods to the solution of problems of high speed 
compressible flow. In this area, the authors [18-20] have proposed a two-step explicit 
implementation of a second order Taylor-Galerkin procedure [16.17] and have used this 
approach to solve successfully a variety of inviscid and viscous problems. The addition 
of artificial viscosity is required to stabilize this solution procedure when it is applied 
to the analysis of problems involving strong discontinuities, and this has the effect of 
spreading flow discontinuities over several computational cells. 

Solution methods based upon high resolution schemes [1-6] give sharper definition 
of flow discontinuities and are supposedly more robust. In two and three dimensions, 
these methods are generally implemented by using operator splitting and applying 

. • • Jti - A 





Here />./>.*. 7 and k denote the density, pressure, specific total energy, temperature 
and thermal conductivity of the fluid respectively and u, is the component of the fluid 
velocity in the direction x, of a Cartesian coordinate system. The equation set is 
completed by the addition of the state equations 

, . 1 ... 1 

P = < '• - * - -MyMyJ • 1 = C,.i# - ->IjUj ; . t:|i 

which are valid for a perfect gas. where *• is the ratio of the specific heats and c,. is the 
specific heat at constant volume. The components of the viscous stress tensor a tJ are 
given by 


du, iht. duu , 

~ M dx, * d* t ] * X dT k *' i 


(4) 


and it is assumed that A and p are related by 




(5) 


THE FLOW SOLVER: FEM-FCT 


As stated above, high resolution, monotonicity preserving schemes must be devel- 
oped in order to be able to simulate the strong nonlinear discontinuities present in the 
flows under consideration. Although the pertinent literature abounds with high reso- 
lution schemes ll-6j. only Zalesak's generalization [7j of the 1-D FC I schemes of Boris 
and Book [S-10j can be considered a truly multidimensional high resolution scheme. We 
remark here that the use of unstructured grids requires such truly multidimensional 
schemes, as the lack of lines or planes in the mesh inhibits the use of operator splitting. 

Erlebacher [11]. and Parrot and Christie [121 first analyzed FCT schemes in the 
context of finite element methods. We develop their ideas further to include the con- 
sistent mass, which yields high temporal accuracy, and to systems of equations. 


The Concept of Flux- Corrected Transport (FC'Ti 

We consider a set of conservation laws given by a system of partial differential equa- 
tions of the form given in eqn.(l). and assume that the advective fluxes F a = F a {1. ) 
play a dominant role over the viscous fluxes F v = F V {F). For flows described by 
eqn.(l), discontinuities in the variables may arise (e.g. shocks or contact discontinu- 
ities). Any numerical scheme of order higher than one will produce overshoots or ripples 
at such discontinuities (the so-called 'Godunov theorem' [15]). Aery often, particularly 
for mildly nonlinear systems, these overshoots can be tolerated. However, for the class 
of problems studied here, overshoots will eventually lead to numerical instability, and 
wiil therefore have to be suppressed. 
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The idea behind FCT is to combine a high-order scheme with a low-order scheme 
in such a way that in regions where the variables under consideration vary smoothly 
(so that a Taylor expansion makes sense) the high-order scheme is employed, whereas 
in those regions where the variables vary abruptly the schemes are combined, in a 
conservative manner, in an attempt to ensure a monotonic solution. 

The temporal discretization ofeqn.(l) yields 

r n+1 = (’" + At r . (6) 

where AI is the increment of the unknowns obtained for a given scheme at time 
t — t m . Our aim is to obtain a A V of as high an order as possible without introducing 
overshoots. To this end. we re-write eqn.(6) as: 

r »+i = V n + AU t + ( {7) 

or 

U m + 1 = V‘ + (A U k - A U‘). (8) 

Here A U h and A U l denote the increments obtained by some high- and low-order scheme 
respectively, whereas U 1 is the monotone, ripple-free solution at time t = t* +l of the 
low-order scheme. The idea behind FCT is to limit the second term on the right-hand 
side of eqn.(8): 


C B+1 = U l -I- /im( AC* - A C'). (9) 

in such a way that no new over /undershoots are created. 

It is at this point that a further constraint, given by the conservation law (1) 
itself must be taken into account: strict conservation on the discrete level should be 
maintained. The simplest way to guarantee this for node-centered schemes (and we 
will only consider those here) is by constructing schemes for which the sum of the 
contributions of each individual element (cell) to its surrounding nodes vanishes ("all 
that comes in goes out'). This means that the limiting process (eqn.(9)) will have to 
be carried out in the elements (cells). 

Algorithmic Implementation 

We can now define FCT in a quantitative way. We follow Zalesak's exposition [7J. 
but modify the term ‘flux’ by ‘element contribution to a node’. Those more familiar 
with finite volume or finite difference schemes should replace ‘element’ by ‘cell’ in what 
follows. 

FCT consists of the following six algorithmic steps: 

1) Compute LEC: the ‘low-order element contribution’ from some low-order scheme 
guaranteed to give monotonic results for the problem at hand: 
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2) Compute HEC: the ‘high-order element 
scheme; 


contribution’. 


given by some high-order 


3) Define AEC: the ‘antidiffusive element contributions’ : 

AEC = HEC - LEC 


4) Compute the updated low-order solution : 

U l = U m + Y. LEC =(/' + AU' 

el 


( 10 ) 


5) Limit or ‘correct’ the AEC so that U”+' 
extrema not also found in or U 9 : 


as computed in step 6 below is free of 


AEC* = Cel * AEC , 0 <Cel< 1; 

6) Apply the limited AEC : 

= U 1 + £ AEC*. 

el 


( 11 ) 

( 12 ) 


The Limiting Procedure 


We ^ ° n S “P 5 *»ove. 

a) Pf : the sum of all positive (negative) antidiffusive element contributions to node I 

''=£{ZW,, 


b) Q r . the maximum (minimum) increment (decrement) node I is allowed 
in step o above 


to achieve 


Qf = Up" - l rl 


"I*— th <“ maximum (minimum) value the un- 
known U at node I is allowed to achieve in step 6 above. 

c) 


R ± ;= / min{\,Q±/p±) if /»+ > o, p~ < 


\ 


0 if P± = o 


Now take, for each element: 
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(13) 


Cel = minUUment nodes) { £ '' ^ > ° 0 


Finally, we obtain If"' in three steps : 

a) maximum (minimum) nodal U of U n and V 1 : 

b) maximum (minimum) nodal value of element : 

{ ' ; ' = {Z} (l - r » **> • 

where A. B C represent the nodes of element el. 

c) maximum (minimum) V of all elements surrounding node I 


rj’z: _ J max ltr- r* 

1 ‘ -{min!' " Lm) ■ 


where 1.2, .... m represent the elements surrounding node I. 

This completes the description of the limiting procedure. Up to this point we have 
been completely general in our description, so that eqns.(6)-( 13) may be applied to 
any FEM-FCT scheme. In what follows, we restrict the exposition to the finite element 
schemes employed in the present work, describing the high and low-order schemes used. 

The High-Order Scheme: Consistent-Mass Tavlor Galerkin 


As the high-order scheme, we employ a two-step form [18-20] of the one-step Taylor- 
Galerkin schemes described in [16,17]. These schemes belong to the Lax-Wendroff 
class, and could be substituted by any other high-order scheme which appears more 
convenient, including implicit schemes. Given the system of equations (1), we advance 
the solution from t n to t n+1 = t n 4- At as follows: 


a) First step (advective predictor): 


(f n+ 2 = U n 


At 

2 


dFJ_ 

dij 


(14) 


b) Second step : 


A U* = U n+1 - U n = -At 


dFf 

dij 


+ At 


dFj 

dx } 


(15) 
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The spatial discretization of (14) and (15) is performed via the classic Galerkin 
weighted residual method (18-20], using linear elements, i.e. 3-noded triangles in 2-D 
and 4-noded tetrahedra in 3-D. For (15) the following system of equations is obtained: 

M c AU n = R n . (16) 

where Me denotes the consistent mass matrix [18-20], AC the vector of nodal incre- 
ments and R the vector of added element contributions to the nodes. As Me possesses 
an excellent condition number, eqn.(16) is never solved directly, but iteratively, requir- 
ing typically three passes (17). We recast the converged solution of eqn.(16) into the 
following form, which will be of use later on : 


M l • AU k * R + {M l - M c ) • AU h . (17) 

Here Ml denotes the diagonal, lumped mass-matrix (see [17]). 

The Low-Order Scheme: Lumped-Mass Taylor Galerkin plus Diffusion 

The requirement placed on the low-order scheme in any FCT-method is mono- 
tonicity. The low-order scheme must not produce any artificial, or numerical, ‘ripples' 
or ‘wiggles’. It is clear that the better the low-order scheme, the easier the resulting 
task of limiting will be. Therefore an obvious candidate for the low-order scheme is 
Godunov’s method [15]. However, this scheme would be relatively expensive, and its 
extension to unstructured grids remains unclear. 

We have so far added ‘mass-diffusion’ to the lumped-mass Taylor-Galerkin scheme 
in the context of FEM-FCT [13,14]. This simplest and least expensive form of diffusion 
is obtained by subtracting the lumped mass-matrix from the consistent mass-matrix 
for linear elements: 


DIFF = C4 • (Me - Ml) • U*. 

The element matrix thus obtained for 2-D triangles is of the form 

{ 2 -i -r 
-1 2 -1 • 
- 1-1 2 


(18) 


(19) 


Observe that we cannot simply add this diffusion to the high-order scheme in order 
to obtain monotonic results, as a multipoint-coupling of the right-hand side occurs due 
to the consistent mass-matrix employed in the high-order scheme . The imposition of 
monotonicity can nevertheless be achieved by using a lumped mass-matrix instead. As 
the terms originating from the discretization of the fluxes F * in ( 1 ) are the same as in 
(15), the low-order scheme is given by 

M l AU‘ = R + DIFF. (20) 



i 
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Resulting Antidiffusive Element Contributions 


Subtracting (20) from (17) yields the equation 

Ml ■ {AV h - AC 1 ) = R + {M l - M c ) ■ AU h - R- D1FF . (21) 

or. using eqn.( 18) 


AU h - AU l = Mr 1 • ( M l - M c ) ■ (e 4 ■ C n + A U h ). (22) 

Note that all terms arising from the discretization of the fluxes F' in ( 1),( 15), (20) 
have now disappeared. This is of particular importance if the antidiffusive element 
contributions must be recomputed (small core memory machines), and real eas effects 

are taken into account (table look-up times are considerable) or real viscosity effects 
have to be included (Navier- Stokes equations). 

Limiting for Systems of Equations 

The results available in the literature (8-10] and our own experience [13.14] have 
shown that with FCT results of excellent quality can be obtained for a single PDE. 
However, when trying to extend the limiting process to systems of PDEs. no imme- 
diately obvious or natural limiting procedure becomes apparent. Obviously, for 1-D 
problems one could advect each simple wave system separately, and then assemble the 
solution at the new time step. However, for multidimensional problems such a splitting 
is not possible, as the acoustic waves are circular. FDM-FCT-codes used for production 
runs (21,22] have so far limited each equation separately, invoking operator-splitting 
arguments. This approach does not always give very good results, as may be seen from 
Sod’s comparison of schemes for the Riemann problem [23], and has been a point of 
continuing criticism by those who prefer to use the more costly Riemann-solver-based, 
essentially one-dimensional TVD-schemes [1-6]. It would therefore appear as attrac- 
tive to introduce ‘system character' for the limiter by combining the limiters for all 


some undershoots for very strong shocks are present. This option is currently our 
preferred choice for transient problems. 

iii) Use of the minimum of the limiters obtained for the density and the pressure 
(C'i = min(C e i(density).C'.i(prest>ure})) : this again produces acceptable results, 
particularly for steady-state problems. 


NUMERICAL EXAMPLES 


a) Shock over an indentation: The first problem considered simulates the transient 
ftowfieid produced by the interaction of a strong shock with an indentation in the 
ground. For this case, the shock Mach number was set to M, = 25, which corresponds 
to a pressure-jump ratio of about 1:100. During the transient, pressure ratios as high 
as 1:1000 result. The problem statement, solution domain, spatial discretization and 
solutions obtained are shown in Figs, la- le. Note that an adaptive refinement scheme for 
transient problems [26) was used to reduce the overall storage and CPU requirements. 

As the shock travels over the indentation, it produces a bow shock and a rar- 
efaction (Figs.la.lb). Then, it collides with the right wall of the indentation and 
bounces back, producing several shock/shock and shock/contact discontinuity interac- 
tions (Figs.lc.ld). Observe the level of physically relevant detail that the scheme is 
able to reproduce, e.g. the triple shock produced at T=0.12 (Figs. Id, le). The veloc- 
ity pattern generated by these interactions has been magnified in Fig.le, and shows a 
large residual vortex, as well as the different shock fronts and other discontinuities. We 
remark that at all times the shocks are captured within 2 to 3 elements. 

In the present case, we used as limiter for all equations the minimum of the limiters 
computed for the continuity and energy equations. It is found, that for the strong shocks 
present in such flowfields, even a pressure-undershoot of 0.1% will lead to negative 
pressures. Therefore, the pressure is additionally limited artificially in order to be 
positive (albeit small) at all times. 

b) Steady supersonic flow past a circular cylinder: This problem involves inviscid 

Mach 3 flow past a circular cylinder. The solution has been obtained by relaxing, 
with local timesteps, the transient solution towards the final steady-state. During this 
iteration process, the grid was adapted three times to the solution by using an adap- 
tive mesh regeneration technique [27] . The final grid is shown in Fig.2a. A detail of 
the pressure coefficient distribution is shown in Fig.2b, and the variation of pressure 
coefficient along the centre line and over the cylinder surface is given in Fig.2c. 

cl Shock-bubble interaction: This problem is included here to demonstrate a new ant- 
isymmetric capability, and also to show that not only geometrically complex domains, 
but also physically complex problems can be handled economically by the methodolo- 
gies developed. Initially, a weak shock (M a = 1.29), coming from the left in Fig.3a. 
travels into a bubble of heavier material. In the present case, the outer medium was 
assumed to be air, while the bubble was assumed to consist of freon. Due to the higher 
density of freon, the shock speed inside the bubble decreases (Fig.3b). While the outer 
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shock bends over, the inner shock focuses at the right end of the bubble, producing a 
significant overpressure (Fig.3c). and initiating a small, circular blast wave (Fig. 3d). 

d ) Steady supersonic flow over a flat plate: The fourth problem considered is the 

steady state solution of supersonic viscous flow over a flat plate. The flow condi- 
tions correspond identically to one of the cases considered by Carter [28], using a finite 
difference scheme. The free stream Mach number is 3 and the Reynolds number based 
on the plate length is 1000. The temperature of the plate is assumed constant. The 
Sutherland viscosity law (see. e.g. Schlichting [29]) is used and the initial conditions 
are chosen to be appropriate to the case of a flat plate impulsively inserted into the free 
stream. The mesh used is displayed in Fig.4a, and the general features of the solution 
can be appreciated in the density contour plots shown in Fig.4b. The variation of the 
computed wall pressure distribution is given in Fig.4c. 

CONCLUSIONS 

It has been demonstrated how unstructured grids and high resolution schemes 
may be combined, yielding FEM-FCT. The numerical examples indicate that a high 
accuracy can be obtained economically for problems involving complex domains and/or 
adaptive mesh refinement. Furthermore, the ‘equation-splitting' employed in classic 
FCT-codes [21,22] has been extended by coupling or ‘synchronizing' the limiters of all 
the equations involved, without taking recourse to more costly Riemann- solver- based 
monotone schemes. 

Extensions of the present work are under investigation and involve the development 
of better limiters for systems of equations in the context of FEM-FCT, the extension of 
FEM-FCT to implicit or semi-implicit time-stepping schemes [31], and the combination 
of FEM-FCT with unstructured multigrid methods [32] for the rapid solution of steady 
state problems. 
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